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Abstract 


The wolf spider Pardosa astrigera L. Koch, 1878, an important biological control agent for pests in agriculture, is widely distributed 
in various ecosystems across East Asia. This study used mitochondrial DNA and aimed to provide an in-depth understanding of popu- 
lation genetic structure and evolutionary history throughout the species. Mitochondrial gene sequences from 107 samples of P. astrig- 
era from 25 East Asian populations were used for genetic analyses. Our data revealed an asymmetric phylogeographic distribution in 
two sympatric lineages (1—2) of P. astrigera in continental East Asia. The spatio-temporal pattern of two mitotypes of P. astrigera in 
this region gives strong support for a Northeast Asian origin during the late Pleistocene (~1.69 million years ago) and the population 
expansion time of ~74,340 (58,832—104,236) years ago (during the last glacial period) and dual colonization around East Asia from 
two directions: from North to South and from East to West. Our phylogeographic results suggested that Pleistocene climate oscilla- 
tions with subsequent fragmentation events and secondary contacts were the major impact factors of the diversification, geographic 
distribution, and expansion patterns of P. astrigera, and human activities and ballooning probably accelerated its recent dispersal. 
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Introduction 


Pardosa astrigera (Araneae: Lycosidae), a wandering 
spider, inhabits a relatively arid-cold environment. It is a 
common species in various ecosystems of the temperate 
and subtropical regions of East Asia and possesses strong 
environmental adaptability, high dispersal potentiality, 
and rapid population diversification (Lida et al. 2016; Li 
et al. 2020). This spider can prey on a variety of pests in 
agriculture and forestry and has strong hunger tolerance 
and high fertility, as well as long-term population stabil- 
ity and substantially large populations (Li et al. 2020). 
Therefore, P. astrigera plays a vital role in pest control. 
It is widespread in East Asian regions such as China, 
Korea, Japan, and Russia (Far East) (Song et al. 1999; 
World Spider Catalog 2024). This wolf spider is a dom- 
inant species in various ecosystems across Korea (Kim 


and Yoo 1997) and most parts of China, except Fujiang, 
Guangdong, and Hainan provinces (Song et al. 1999). 
Previous morphological studies suggested that P. astrig- 
era was likely a species complex (Schenkel 1963). Re- 
cent investigations have considered the morphological 
variations among the different geographical populations 
to be intraspecific (Yin et al. 1997; Chang et al. 2007). 
Studies on the population genetic structure of P. astrigera 
can estimate genetic diversity and gene flow among its 
different populations and infer the ancestral populations 
and the phylogeographic pattern. Such information is es- 
sential for effectively protecting and utilizing the natural 
predators of pests in agriculture and forestry. 

The phylogeography of P. astrigera distributed in lo- 
cal regions of China was investigated using nuclear ITS2 
gene sequences, and two major lineages were observed 
(Chang et al. 2007). However, for the ancient spider 
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species with an extensive geographic range, the infor- 
mation concerning genetic structure, spatio-temporal 
evolution patterns, and demographic history across the 
whole distribution range in East Asia was still unclear, 
and the following questions concerning this widespread 
species have not been investigated: 1) genetic relation- 
ships among the East Asian populations based on mito- 
chondrial genes; 2) allopatric divergence among the East 
Asian regions; and 3) origin and dispersal across East 
Asia. These questions will provide a better understanding 
of the P. astrigera population history in East Asia. 

Because of the haploid maternal inheritance of mito- 
chondrial genes and their high mutation rate and abun- 
dance in cells, they are widely used as molecular markers 
to analyze the population genetic structure and phylo- 
geography of broadly distributed animal species (Avise 
2000). To assess the mtDNA population genetic structure 
and potential dispersal routes of P. astrigera, we studied 
its genetic structure and phylogeography using mitochon- 
drial COI, 16S, and NADH] genes. We recently collected 
individuals of the spider from China and augmented our 
sampling with sequences from GenBank covering nearly 
the entire geographic range of P. astrigera. The current 
data set covers 25 distinct localities (21 from China, three 
from South Korea, and one from Japan) around East Asia. 
It thus offers a more holistic insight into P. astrigera ma- 
triline history than have prior studies. 


Materials and methods 
Sampling, sequencing, and sequence analyses 


We collected 19 samples of P. astrigera from 6 provinces 
of China (Suppl. material 1: table S1) and sequenced 3 
mitochondrial genes (COI, 16S, and NADH7) using three 
primer pairs (Suppl. material 1: table S2) according to the 
procedures of Luo and Li (2022). Additional sequences of 
the species from Korea, Japan and 14 provinces of Chi- 
na were taken from Chang et al. (2007), Li et al. (2011), 
and GenBank (Suppl. material 1: table S1). The mtDNA 
sequences of 107 colonists from 25 geographic sites in 
East Asia were used for genetic analyses (Fig. 1a; Suppl. 
material 1: table S1). BioEdit was used to check the qual- 
ity of the DNA sequences (Hall 1999). Sequence align- 
ments were done with MAFFT v7 (Katoh and Standley 
2013). Unique haplotypes were identified using DnaSP 
v5 (Librado and Rozas 2009). The data matrix is avail- 
able from the authors and will be submitted to the Dryad 
database (online at http://datadryad.org/). 


Phylogeny and haplotype network reconstruction 


A matrix combining the sequences of the COZ, 16S, and 
NADHI markers was constructed and used to obtain the 
phylogenetic relationships of 107 P. astrigera individ- 
uals using the Bayesian inference (BI) and maximum 
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likelihood (ML) approaches. Pardosa laura was used 
as an outgroup. Bayesian analysis was performed in 
MrBayes V3.2.1 (Ronquist and Huelsenbeck 2003). We 
used jModelTest to select the best-fit evolutionary mod- 
el for each codon base or gene partition (Posada 2008) 
with the Bayesian information criterion (Suppl. material 
1: table S3). The Markov chain Monte Carlo (MCMC) 
simulation was run for 15 million generations and sam- 
pled every 1,000 generations, and unlinked parameters 
among partitions were used in the chain runs. Tracer 
v1.5 was employed to check effective sample size (ESS) 
values (>200) for each parameter (Rambaut and Drum- 
mond 2009), and a 50% majority rule consensus tree was 
computed after discarding the first 25% of sampled trees 
as burn-in. The ML tree was reconstructed using IQ- 
Tree v1.6.12 (Nguyen et al. 2015) as implemented in the 
W-IQ-Tree web server (Trifinopoulos et al. 2016). The 
substitution model for each codon base or gene partition 
was automatically selected, and an ultrafast bootstrap 
with 1,000 iterations was conducted to estimate nodal 
support across the topology (Minh et al. 2013; Kalyaana- 
moorthy et al. 2017; Hoang et al. 2018). To compare ge- 
netic connections between geographic populations, three 
haplotype networks were inferred in PopART (Leigh and 
Bryant 2015) using unrooted minimum-spanning trees 
(Bandelt et al. 1999) based on the separate mtDNA gene 
sequences. Based on climate, geographical divisions, 
distributional data for the P. astrigera species from the 
literature (Song et al. 1999; World Spider Catalog 2024), 
and information from our samples. We delimited 6 areas, 
including Northwest China (Qinghait+Gansu, China); 
Northeast Asia (Neimenggu+Heilongjiang, China+Ko- 
reat+Japan); East China (Jiangsu+Shandong+Zheyiang, 
China); Southwest China (Yunnan+Guizhou, China); 
Central China (HubeitHunan+Anhui, China); and North 
China (Hebei+Henan+Shanxit+Sanxi, China). 


Genetic diversity and distance calculation 


We defined six populations from the sampled regions for the 
haplotype network reconstruction above. Population genet- 
ic diversity was inferred from the CO/, 16S, and NADH] 
sequences. Estimations of genetic diversity (number of 
haplotypes, nh; nucleotide diversity, 2; haplotype diversity, 
h) from each population were performed in DnaSP 5.10.01 
(Librado and Rozas 2009). Genetic distances among hap- 
lotypes and overall mean distances of P. astrigera popula- 
tions were estimated using MEGA v5 (Tamura et al. 2011). 


Estimations of divergence time and mtDNA 
substitution rate 


Divergence times and mtDNA substitution rates 
were inferred based on the combined CO/, 16S, 
and NADHI sequences. They were implemented in 
BEAST v1.6.1. (Drummond and Rambaut 2007). 
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Figure 1. a. Sampling regions of the East Asian wolf spider Pardosa astrigera and the haplotype distribution of the /6S gene. 
Detailed sampling information is presented in Suppl. material 1: table S1 in the Supplemental materials; b. Time-tree of the species 
Pardosa astrigera based on the combined mitochondrial CO/, 16S, and NADH] sequences. Blue dots mark the four fossil calibra- 
tion nodes. Numbers at nodes indicate the main divergence times of the haplotypes of the species. “*” indicates stable branches with 
Bayesian support >0.90. 
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First, a test was run to determine if a molecular clock mod- 
el was appropriate based on the harmonic mean (HM) of 
the likelihood values of the MCMC samples (Newton and 
Raftery 1994) and a more accurate stepping-stone (SS) 
sampling method (Xie et al. 2011) in MrBayes v3.2.1. We 
found that the total data sets did not conform to a strict 
molecular clock (P<0.001) and that a relaxed molecular 
clock method was more suitable for estimating divergence 
time. Therefore, we estimated divergence times using the 
uncorrelated lognormal relaxed molecular clock model. 
The tree prior was set to the birth-death and Yule specia- 
tion process models. We repeated the analysis with birth- 
death and Yule prior to assess the sensitivity of our re- 
sults to tree-prior specification. According to the marginal 
likelihood estimated using stepping-stone sampling, the 
birth-death model outperformed the Yule model. Hence, 
subsequent analyses focused on the chronogram result- 
ing from the birth-death prior. }ModelTest (Posada 2008) 
was used to select the best-fitting substitution model for 
each gene or codon base partition under the Bayesian in- 
formation criterion (Suppl. material 1: table S4). To ob- 
tain reliable results, we ran five independent MCMC tree 
searches for 300 million generations and sampled every 
1,000 generations. The convergence of MCMC chains 
was assessed in TRACER v1.5 (Rambaut and Drummond 
2009). TreeAnnotator v1.6.1 (in the BEAST package) 
was used to produce a maximum clade credibility (MCC) 
tree with median heights after discarding 25% of the trees 
as burn-in. We dated the tree of all sampled P. astrigera 
haplotypes from East Asia. For the molecular clock anal- 
ysis, we used the minimum ages based on fossils of Ly- 
cosidae (15 million years ago, Ma; Iturralde-Vinent and 
MacPhee 1996), Oxyopidae (43 Ma; Wunderlich 2004a; 
Magalhaes et al. 2020), Thomisidae (43 Ma; Wunderlich 
2004b; Magalhaes et al. 2020), and Selenopidae (53 Ma; 
Penney 2006; Magalhaes et al. 2020) as calibration points 
(Suppl. material 1: table S4; Renner 2005; Donoghue and 
Benton 2007). The outgroups include the Agelenidae, 
Thomisidae, Oxyopidae, Psechridae, Trechaleidae, Sele- 
nopidae, and some Lycosidae species (Suppl. material 1: 
table S4). Their gene sequences were available from Gen- 
Bank (Suppl. material 1: table S1). 


Demographic analysis 


To test for range expansions, Fu’s FS, Tajima’s D, and 
the mismatch distribution were calculated (Tajima 1989; 
Fu 1997; Jaeger et al. 2005; Smith and Farrell 2005) in 
Arlequin v3.5 (Excoffier and Lischer 2010), with 10,000 
permutations to test for significance. A significant, large 
negative value for Fu’s FS, a significant value for Taji- 
ma’s D, or an unimodal shape of the mismatch distribu- 
tion indicates population expansion (Slatkin and Hudson 
1991; Rogers and Harpending 1992; Aris-Brosou and Ex- 
coffier 1996; Tajima 1996; Fu 1997). 

We used a mean substitution rate (V) to infer the expan- 
sion time of P. astrigera. The value of V was estimated from 
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the above mitochondrial molecular clock analysis. The gen- 
eration time in East Asia for P. astrigera is 6 months (Chen 
et al. 2010; Yang et al. 2018). Therefore, the substitution 
rate per generation was given as v =0.5V. Once v was esti- 
mated, the coalescence time in generations (7) was obtained 
by the relationship t= t/(2 x 2 mv) (Rogers and Harpending 
1992; Harpending et al. 1993; Rogers 1995), where m was 
the concatenated CO/, 16S, and NADH] sequence length, 
and the value of t was calculated using Arlequin v3.5. 

Bayesian skyline plots (BSPs) reconstruct historical 
population sizes from mtDNA genealogies (Drummond 
and Rambaut 2007). BSPs were implemented in BEAST 
v1.6.1. The substitution rate (2.58%) of the combined 
mitochondrial gene sequences was estimated from the 
mitochondrial molecular clock analysis. The substitution 
models were selected using }ModelTest (Suppl. material 
1: table S4). For Bayesian skyline coalescent tree priors, 
we chose the piecewise-linear skyline model. Otherwise, 
default parameters were used. The chain was run for 10 
million steps under an uncorrelated lognormal relaxed 
clock model and sampled every 1,000 steps. The result 
of the Bayesian skyline plot was checked and analyzed 
using Tracer v1.5 with a burn-in of 10%. 

The expansion history of P. astrigera was increasingly 
estimated using Bayesian binary MCMC analyses (BBM, 
a method considered the most general and complex model 
in biogeographical reconstruction; Sanmartin et al. 2001). 
BBM was implemented in RASP v3.0 (Yu et al. 2014) 
using a maximum of two areas per node. The 10,000 trees 
obtained from the BEAST analysis were input into RASP. 
The 1,000 random trees were used for the biogeographi- 
cal reconstruction (BBM). The MCC tree produced in the 
BEAST analysis was used as the input tree. Based on the 
dated chronogram for P. astrigera, we pruned the tree to 
remove outgroups. We divided the range into the same six 
biogeographical areas as above. 


Results 
Sequences 


The 96, 81, and 76 sequences were obtained for the COL, 
16S, and NADH] genes from all samples of P. astrigera, 
respectively. The aligned data of the sequences from P. as- 
trigera had lengths of 932 base pairs (bp) for COI, 574 bp 
for 16S, and 584 bp for NADH. In total, 62, 17, and 38 
haplotypes were found for the COJ, 16S, and NADH7 gene 
sequences, respectively. Among the 62 unique COI haplo- 
types, Hap20 dominated in 8.3% of the samples, Hap39 
in 6.3%, and each of Hap3, Hap15, and Hap52 in 4.2%; 
among the 17 unique /6S haplotypes, Hap4 dominated in 
55.6% of the samples, Hap15 in 14.8%, and Hap9 in 7.4%; 
and among the 38 unique NADH haplotypes, Hap13 dom- 
inated in 36.8% of the samples, Hap22 in 6.6%, and Hap] 
in 5.3%. The concatenated mitochondrial genes comprised 
2,090 bp. All sequences are deposited in GenBank (for ac- 
cession numbers, see Suppl. material 1: table S1). 
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Molecular phylogeny 


In the dated phylogenetic tree, P. astrigera is composed 
of the two lineages (1-2; Fig. 1b). Lineage 1 consists of 
the sample (CF5) from Neimonggu, China, the four sam- 
ples (LY1—3 and LYS) from Henan, China, and the sam- 
ple (MS3) from Heilongjiang province of China, and the 
lineage spiders occur on a clade of the ML tree (Suppl. 
material 1: fig. S1b); and Lineage 2 consists of all other 
individuals of P. astrigera across continental East Asia, 
and the lineage samples cluster together in a clade of the 
BI tree (Suppl. material 1: fig. Sla). The BI and ML trees 
(Suppl. material 1: fig. S1) indicated five sister groups 
(2a—e). The samples in Group 2a are from Neimenggu, 
Shanxi, and Zhejiang provinces of China and South Ko- 
rea; Group 2b from Heilongjiang and Zhejiang provinces 
of China and Japan; Group 2c from Shandong, Gansu, 
and Hubei provinces of China; Group 2d from Sanxi 
and Zhejiang provinces of China; and Group 2e from 
Neimenggu, Hebei, Shandong, and Shanxi provinces of 
China and South Korea. However, the relevance of these 
groups is not clear in the molecular phylogeny. 


Haplotype network 


The haplotype network analyses based on both 16S and 
NADHI genes revealed clear genetic structuring among 
the East Asian populations of P. astrigera (Fig. 2a, b); 
in contrast, those based on the CO/ gene showed a com- 
plex and intricate haplotype connection pattern (Fig. 2c). 
In 16S and NADH] networks, all studied individuals 
clustered into two lineages (1-2; Fig. 2a, b). Lineage 1 
comprises three haplotypes from Neimenggu and Henan 
regions in China; Lineage 2 is composed of the 35 and 14 
haplotypes from all analyzed regions, with one (Hap13) 
and two (Hap4 and Hap15) universal haplotypes occur- 
ring at the center of a star-like cluster of rare haplotypes 
in the NADH/ and 16S networks, respectively (Fig. 2a, 
b). The CO/ network structure consists of 62 haplotypes 


with Hap20 at the center, indicating that the wolf spider 
has complex population genetic relationships and a high 
gene flow among populations (Fig. 2c). 


Genetic diversity 


Northeast Asia and North China populations exhibited 
higher genetic diversity than the other four regional pop- 
ulations based on CO/ and NADH genes (Table 1). The 
genetic diversity of the East China population was higher 
than that of the Central China population based on the 
two protein-coding genes (Table 1). Among all analyzed 
populations, the population of Southwest China had the 
lowest genetic diversity based on CO/ and 16S genes, 
whereas Northwest China had the lowest genetic diversi- 
ty based on the NADH7 gene (Table 1). 


Haplotype divergence times and population 
expansion time 


Based on fossil calibrations (Suppl. material 1: table S5), 
the estimated initial haplotype divergence of P. astrigera 
started at ~1.69 Ma (95% credibility interval, CI: 0.95— 
2.66 Ma; Fig. 1b). Diversification of the five haplotypes 
(Lineage 1) from Heilongjiang, Neimenggu, and Henan 
provinces of China occurred ~0.98 Ma (95% CI: 0.50-1.44 
Ma; Fig. 1b). The two haplotypes from Jiangsu and Zheji- 
ang provinces of China diverged from the other haplotypes 
of Lineage 2 at ~1.21 Ma(95% CI: 0.62-1.95 Ma; Fig. 1b). 
The mean substitution rate (V) of the concatenated mtDNA 
sequences for the species was estimated to be 0.0129 (95% 
CI: 0.0092—0.0163) per million years, equivalent to a di- 
vergence rate of 2.58% (1.84%-—3.26%) per million years. 
The t value was inferred from the mitochondrial allele 
frequencies and linkage disequilibrium to be 2.244. Using 
the values of V and t, the estimated generation time since 
expansion for the P. astrigera populations in East Asia was 
approximately 74,340 (58,832—104,236) years ago. 


Table 1. Summary of genetic diversity in P. astrigera from East Asia based on mtDNA sequences. The number of haplotypes (nh), 
nucleotide diversity (2), haplotype diversity (h), Fu’s FS, and Tajima’s D are shown. Northeast Asia (Neimenggu+Heilongjiang, Chi- 
na+Korea+Japan); Northwest China (Qinghait+Gansu, China); North China (Hebeit+Henan+Shanxi+Sanxi, China); Central China 
(Hubeit+Hunan+Anhui, China); East China (Jiangsu+Shandong+Zhejiang, China); and Southwest China (Yunnan+Guizhou, China). 


Population "NADH... 
i Fu’s FS Tajima’sD ni nh TT 


Col 
ni nh TT h 
Northeast 10 9 0.014+ 0.978+ -3.323 -0.90274** 9 4 0.004 + 
Asia 0.003 0.054 0.001 
Northwest 11 10 0.007+ 0.982+ -6.904 -1.69129** 10 6 0.005+ 
China 0.002 0.046 0.001 
North China 27 23 0.010+ 0.986+ -21.022 -1.39856** 19 7 0.003 + 
0.001 0.015 0.001 
Central 13 10 0.009+ 0.923+ 4.517 -1.20786** 10 6 0.004 + 
China 0.002 0.069 0.001 
East China 23 21 O.011+ 0.992+ -18.844 -1.51301** 22 7 0.003 + 
0.001 0.015 0.001 
Southwest 12 7 0.006+ 0.879+ -1.828 -1.25306** 11 3 0.002 + 
China 0.002 0.075 0.001 
Total 96 58 0.010+ 0.979 + -72.764 -1.58133** 81 17 0.003 + 
0.001 0.006 0.000 


16S NADH1 

h Fu’s FS Tajima’sD ni nh TT h Fu’s FS_ Tajima’s D 
0.694+ 0.722 -0.84257** 9 8 0.010+ 0.972+ 4.093 -0.95811** 
0.147 0.002 0.064 
0.889 + -2.781 -0.78138** 10 6 0.005+ 0.778+ -2.521 -1.57285** 
0.075 0.001 0.137 
0.544+ 4.179 -1.20300** 19 11 0.009+ 0.895+ -3.876 -1.45647** 
0.136 0.002 0.057 
0.844+ -3.412 -1.38818** 9 6 0.005+ 0.833+ -2.495 -1.79752* 
0.103 0.002 0.127 
0.671 + -3.870 -1.46068** 18 11 0.005+ 0.856+ -7.230 -2.05890* 
0.094 0.001 0.079 
0.473 + 0.659 -0.77815** 10 7 0.006+ 0.867+ -3.347 -1.68719** 
0.162 0.002 0.107 
0.668 + -15.586 -1.85939* 76 38 0.007+ 0.860+ -45.933 -2.39099* 
0.054 0.001 0.039 


(East Asia) 
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Figure 2. a. /6S haplotype network of the Pardosa astrigera spiders; b. Network from NADH]; ¢c. Network from COZ: NHJ, 
Northeast Asia (Neimenggu+Heilongjiang, Chinat+Korea+Japan); SGS, North China (Hebei+Henan+Shanxi+Sanxi, China); QHS, 
Northwest China (Qinghait+Gansu, China); YUG, Southwest China (Yunnan+Guizhou, China); HUA, Central China (Hubei+Hu- 
nan+Anhui, China); JSZ, East China (Jiangsu+Shandong+Zheyjiang, China). 
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Demographic history 


The neutrality test of mitochondrial genes showed all six 
separate populations, and the total population had nega- 
tive values of Fu’s FS/Tajima’s D or significantly nega- 
tive values of Tajima’s D (Table 1). These values indicate 
historical population size expansions, genetic hitchhik- 
ing, and/or selection (Tajima 1989; Fu 1997). The 16S 
and NADH] network structures in which one or two 
universal haplotypes occurred at the center of a star-like 
cluster of rare haplotypes indicated rapid demographic 
expansion of the P. astrigera population in East Asia (Fig. 
2a, b). The unimodal mismatch distribution and its strong 
bias toward low divergence values, as distinguished by 
O- and 1-nucleotide changes, illustrated relatively recent 
expansion from a small number of ancestors (Rogers and 
Harpending 1992) or a range expansion with high levels 
of migration between neighboring demes (Ray et al. 2003; 
Excoffier 2004; Fig. 3a—c). The BSP analyses revealed 
that the East Asian P. astrigera adequate population size 
underwent four stages: 1) a gradual growth since ~50,000 
years ago; 2) a slight fall between 8,000 and 1,800 years 
ago; 3) asharp decrease from 1,800 to 600 years ago; and 
4) arapid rise after ~600 years ago (Fig. 3d). 

The possible ancestral ranges and dispersal pathways 
of P. astrigera in East Asia were inferred by the BBM 
analyses (Fig. 4a). The most likely ancestral area of this 
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spider was Northeast Asia, and it possibly dispersed 
from North China and East China to Northwest China 
or Southwest China via Central China (Fig. 4b). These 
results were further supported by the genetic diversity 
estimations (Table 1) because ancestral populations pos- 
sess higher genetic diversity than derived populations 
(Savolainen et al. 2002). 


Discussion 


Genetic variations of P. astrigera 


Previous studies suggested that P. astrigera was a species 
complex that showed high intraspecific and interspecif- 
ic morphological variations (Schenkel 1963; Yin et al. 
1997; World Spider Catalog, 2024). Our analyses of the 
mitochondrial CO/, 16S, and NADH] data also found the 
largest genetic distance (p-distance) among the P. astrig- 
era haplotypes was 0.033 for COI, 0.013 for 16S, and 
0.030 for NADHT1], and the largest overall mean distance 
(p-distance) among the six populations was 0.014 for 
COT, 0.005 for /6S, and 0.010 for NADH] (Suppl. ma- 
terial 1: table S6). Our DNA data indicated that the wolf 
spider has a rich genetic diversity, and the six divided geo- 
graphical populations showed obvious genetic variations 
(Table 1). Moreover, in the study of Chang et al. (2007), 
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Figure 3. Mismatch distributions of Pardosa astrigera from the whole sample based on the CO/ sequences (a) and based on the 
16S sequences (b) and the NADH sequences (c) independently. Bayesian skyline plots from the whole sample (d); middle lines 
represent median estimates of the effective population size, and shaded areas represent 95% of the highest posterior densities (95% 
HPD). The effective population size is presented on a logarithmic scale. The LGM represents the last glacial maximum. 
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Figure 4. a. Biogeographical reconstruction of P. astrigera; b. Probable dispersion routes in East Asia. Northeast Asia (Neimenggu+Hei- 
longjiang, China+Korea+Japan); North China (Hebei+Henan+Shanxi+Sanxi, China); Northwest China (Qinghait+Gansu, China); South- 
west China (Yunnan+Guizhou, China); Central China (Hubeit+Hunan+Anhui, China); East China (Jiangsu+Shandong+Zheyiang, China). 


two phenotypes (types A and B distinguished by the 
shape and length of the median apophysis in palpus) of 
the P. astrigera males were recognized in China, and they 
differed genetically in the presence (type A corresponding 
to Lineage 1 of this study) or absence (type B correspond- 
ing to Lineage 2) of common insertions and deletions in 
the nuclear ITS2 gene. Therefore, our mtDNA data also 
supported the two differentiated phenotypes. 
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The evolutionary potential of species under envi- 
ronmental stress depends on levels of genetic diversity 
(Black et al. 1992; Gurdebeke et al. 2003). High genetic 
diversity is an important factor in adaptation to environ- 
mental changes. The wolf spider, P. astrigera, is wide- 
spread in East Asia. As a major natural predator of pests 
in various semi-arid land ecosystems such as farmland, 
forests, vegetable fields, tea gardens, and fruit orchards, 
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it has a wide distribution, a large population size, strong 
adaptability, and a complex mechanism for protecting 
offspring (spiderlings clustering on the abdomen and 
back of mothers after hatching and dispersing after one 
molt). Therefore, this wolf spider seems to be highly re- 
sistant to environmental changes and has the genetic basis 
to become a dominant and widespread species. 


Haplotype relationships and mitochondrial genes 


The complex haplotype relationships among popula- 
tions of P. astrigera and the high frequency of gene flow 
within regions were revealed by the CO/ haplotype net- 
work analyses. However, the inferred 16S and NADHI 
networks show clear haplotype relationships and a lower 
degree of divergence than that obtained using CO/. These 
findings can be explained by the importance of CO/ in 
cells and strong natural selection (Arnold 2012). The hap- 
lotype network structure revealed by the mitochondrial 
16S or NADH] analyses is consistent with that by nuclear 
ITS2 analyses (Chang et al. 2007), and we found a similar 
mutation rate of 16S to ITS2 in P. astrigera. 

Both haplotype networks inferred from 16S and 
NADH] are composed of two lineages (1-2). Combined 
with the results of the biogeographical reconstruction 
(BBM), we speculated that the Lineage 1 spiders were 
the ancestors of P. astrigera. The main haplotypes (Hap 4 
and Hap 15 for 16S; Hap 13 for NADH7) at the center of 
Lineage 2 of each haplotype network were identified in 
all populations. Therefore, we proposed them as the an- 
cestral haplotypes within the lineage. Further, the number 
of mutation steps suggested that the other haplotypes of 
Lineage 2 might have derived from the main haplotypes 
over different time periods. 


Phylogeographical pattern and historical 
demography 


Our mtDNA data support that P. astrigera comprises two 
sympatric lineages (1-2). Lineage 1 has a limited geo- 
graphic distribution north of the Qinling Mountains (Mts) 
of China, including the sample (CF5) from Neimonggu 
province of China, the four samples (LY1-3 and LY5) 
from Henan of China, and the sample (MS3) from Hei- 
longjiang of China; in contrast, Lineage 2 is present across 
the entire ranges of the species in continental East Asia. 
This population structure pattern supports the hypothesis 
that intraspecific phylogeography involves common Iin- 
eages that are widespread, plus related lineages that are 
confined to one or a few nearby locales (Avise 2000). The 
two lineages of P. astrigera early split ~1.69 (0.95—2.66) 
Ma, diversified since the late Pleistocene in East Asia, 
and exist sympatrically at multiple locations. The results 
suggest that the Pleistocene glaciation is probably a major 
driving factor in its genetic differentiation. One haplotype 
(Hap 13) of NADH] and two (Hap 4 and Hap 15 of 16S) 
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occurred in all six sampled populations across continental 
East Asia, indicating recent frequent dispersal of the wolf 
spider, probably with human activities. Pardosa astrig- 
era 1S a dominant species and has strong environmental 
adaptability in various open ecosystems such as farmland, 
grassland, urban regions, and hills with human activity. 
Moreover, continental East Asia has a dense population. 
Therefore, human movements and trades can change the 
geographic distribution range of P. astrigera and acceler- 
ate its dispersal across the East Asian region. Additional- 
ly, ballooning is likely another driver for the frequent and/ 
or long-distant dispersal of P. astrigera across East Asia. 
Dispersal by ballooning on silken threads is a well-known 
behavior in wolf spiders (Matthew 1982; Postiglioni et 
ale 2010). 

Analyses of the mtDNA population genetics revealed 
that P. astrigera may have undergone a recent demograph- 
ic expansion in East Asia. Our estimation of the popu- 
lation expansion time [~74,340 (58,832—104,236) years 
ago] of the spider based on fossil calibrations is during 
the last glaciation. The time is earlier than that (28,000 
years ago) inferred by Chang et al. (2007) using a cali- 
bration assuming a rate of 2.3% sequence divergence per 
million years for arthropod mtDNA. Our estimated diver- 
gence rate (2.58%) was comparable to those estimated by 
Papadopoulou et al. (2010), who, through an extensive 
survey of tenebrionid beetles, obtained a divergence rate 
of 2.69% per million years for their concatenated CO/ 
and 16S sequences. 

Lycosidae spiders occur globally, with many species 
reported to have widespread distributions and inhabit 
various ecosystems (Song et al. 1999; World Spider Cat- 
alog, 2024). The wolf spider P. astrigera lives mainly in 
the East Asian regions (Song et al. 1999; World Spider 
Catalog 2024; Fig. la; Suppl. material 1: table S1). Our 
analyses provide novel insights into the diversification of 
P. astrigera in continental East Asia during the Cenozoic. 
This spider is thought to originate most likely in North- 
east Asia around the Pliocene and possibly dispersed in 
the major directions from North to South and from East 
to West. Our findings indicated that the pattern of genetic 
diversity within populations seems predominantly deter- 
mined by historical factors but is modified by contempo- 
rary aspects. 


Impacts of the Pleistocene glaciation on the 
evolution of P. astrigera 


Many species experience retreats, expansions, and di- 
versification in response to cold-warm climatic oscil- 
lations during the Pleistocene gracial periods (Hewitt 
2000, 2004; Ding et al. 2011). They often retracted to 
their southern range or refugia during the cold glaciation 
periods and expanded their distribution ranges during 
the warm interglacial periods, and thus their glacial re- 
fugia are often located at the south and postglacial ex- 
pansions are often northward (Fu and Wen 2023). These 
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evolutionary patterns resulted in a genetic diversity pat- 
tern of “southern richness to northern purity” (Hewitt 
2000). That was, however, not supported by the present 
study, which found the Northeast Asian P. astrigera pop- 
ulation had the highest genetic diversity among the six 
divided geographic populations from continental East 
Asia. Moreover, the ancestral area reconstruction of this 
Species suggests its ancestors inhabited Northeast Asia 
and subsequently dispersed southwestward. Therefore, 
the Northeast Asian regions were probably the glacial 
refugia of P. astrigera. The refugia affected the phyloge- 
ography and genetic structure of the species. That is sim- 
ilar to the results of Gonzalez-Trujillo et al. (2016), who 
studied the Pleistocene refugia of Pardosa sierra on the 
Baja California Peninsula. On continental East Asia, the 
Pleistocene glacial-interglacial climatic cycle plays a cru- 
cial role in some species (Fu and Wen 2023). Pleistocene 
climatic fluctuations have also been shown to influence 
the distribution patterns of European and North Ameri- 
can lycosids (Muster and Berendonk 2006; Ivanov et al. 
2023). These results and our findings suggested that simi- 
lar interglacial dynamics shaped the modern boreal fauna 
of the wolf spiders in the Holarctic region, and secondary 
contacts were likely to be common among the allopatric 
cold-adapted species. 

The spiders of P. astrigera started during the mid-Pleis- 
tocene, resulting in Lineage | and Lineage 2, and subse- 
quently both the two lineages rapidly diversified. There- 
fore, the mid-Pleistocene glaciation is an important driver 
for the diversification of this spider. The population ex- 
pansion of P. astrigera began during the last glaciation. 
We speculated that the co-existence of Lineage 1 and 
Lineage 2 occurred before the last glacial epoch. The last 
glaciation triggered the expansion of the Lineage 2 pop- 
ulations into the south of the Qinling Mts, whereas the 
Lineage | populations were likely blocked by the Qinling 
Mts. Thus, multiple glaciations during the Pleistocene 
affected the arid-cold distribution and diversification of 
P. astrigera in East Asia. 


Future directions 


This study has allowed understanding of the matrilin- 
eage structure of the wolf spider P. astrigera and gener- 
ated hypotheses regarding its origin and dispersal using 
the mitochondrial CO/, 16S, and NADH loci. However, 
while mtDNA data have proven highly useful in phylo- 
geographic analyses, there are several shortcomings to 
this approach (e.g., Avise 2000). There is a potential mis- 
match between the dispersal histories of males and fe- 
males and a likely mismatch between a species tree and 
any single gene tree. In this paper, we summarize avail- 
able COI, 16S, and NADH data, add new data analyses 
of patterns across East Asia, and discuss hypotheses on 
the dispersal history of P. astrigera. Robustly testing the 
historical expansion hypotheses and gaining a deeper un- 
derstanding of the population structure of this spider will 
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require DNA data from the whole genome, using methods 
such as RAD-seq. Such data will tease apart male and 
female dispersal, testing the “female-based” dispersal 
hypotheses proposed here. They will also allow testing 
of the “independent lineages” hypothesis through precise 
measures of genetic isolation and gene flow unavailable 
to several marker studies. Sampling is the process of 
choosing a subset of a target population that will serve as 
its representative. In this paper, the total sampling speci- 
mens for some regions are on the low side. However, we 
hope our study can aid in strategic resampling, reflecting 
known lineage divergences. 
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